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Abstract 

A method for generating angular forces around cr-bonded transition metal 
ions is generalized to treat 7r-bonded configurations. The theoretical approach 
is based on an analysis of a ligand-field Hamiltonian based on the moments of 
the electron state distribution. The functional forms that are obtained involve 
a modification of the usual expression of the binding energy as a sum of ligand- 
ligand interactions, which however requires very little increased in CPU time. 
The angular interactions have simple forms involving sin and cos functions, 
whose relative weights depend on whether the ligands are a- or 7r-bonded. 
They describe the ligand-field stabilization energy to an accuracy of about 
10%. The resulting force field is used to model the structure of small clusters, 
including fragments of the copper blue protein structure. Large deviations 
from the typical square copper coordination are found when 7r-bonded ligands 
are present. 

Running title: Angular Forces around Ions 

Keywords: Metalloproteins, Enzyme structure, Ligand-field splitting, 
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I. INTRODUCTION 



Atomic-level simulations of proteins interacting with transition-metal ions have the 
promise of elucidating a wide variety of biophysical phenomena (Kaim, 1994). Such simu- 
lations can help establish both the native structure and reaction paths of metalloproteins, 
and explain the roles of metal ions in protein folding and stability. Similar applications can 
be seen for simulations of metal ions interacting with DNA and RNA. However, the utility 
of such simulations depends critically on the availability of accurate but computationally 
tractable force fields for metals interacting with proteins. There are no quantitatively accu- 
rate force fields for these purposes, and it is probably not possible to construct force fields 
of a simple enough form to be useful that have quantitative accuracy. However, one can 
hope to generate force fields that include the basic physical effects at a level of accuracy 
sufficient to ascertain chemical trends as metals or ligands vary. The situation in this regard 
is best for "simple" metal ions that have no partly occupied shells. Here an ionic picture 
based on electrostatics, supplemented by empirical repulsive forces, may hope to treat the 
most important physical effects. Even here though, care should be taken, since the bond 
is not completely ionic but has some partial covalent character. For transition metals, the 
situation is much more difficult. The ligand- field splitting of the ci-shell leads to important 
electronic contributions which cannot be ignored. These are manifested, for example, in 
the typical square or tetragonal coordination of Cu 2+ complexes. Such effects cannot be 
described by radial interactions, but rather require the introduction of angular forces de- 
scribing the energetics of the transition-metal ci-shell. Most existing methods for including 
angular forces in simulations of metal ions have used assumed functional forms (Comba et 
al., 1995, Allured et al., 1991, Timofeeva et al., 1995, Wiesemann et al., 1994, Sayle et al., 
1995, Sayle et al., 1997) based, for example on the observed structures of small complexes. 
A recently developed method (Landis et al., 1998) treats the environment-dependent ener- 
getics of transition metals via a valence-bond approach. This method appears promising 
for covalently bonded metals. However, the bonding state of transition metals in proteins 
is may be rather different from the sd n hybridized configuration assumed in that work. In 
the present analysis, we seek to develop a method suitable for more ionic configurations. 

We have recently shown (Carlsson, 1998) how the ligand-field splitting effects of a-bonded 
transition metal ions can be described by a force field of a fairly simple form. One starts 
with an explicitly quantum-mechanical form for the electronic energy, and then then solves 
the quantum-mechanical problem with systematic approximation methods to extract the 
real-space description of the electronic bonding energy. A second-order treatment of the 
hybridization terms between the ligand and metal orbitals is used to generate a ligand- 
field Hamiltonian for the ci-electrons. The electronic bonding energy of this Hamiltonian is 
analyzed in terms its moments. The moment analysis of the energy gives a "semiclassical" 
energy function with a simple trigonometric angular form. It is not precisely an additive 
function of ligand-ligand interactions, but is nearly as simple computationally. This energy 
function gives quite accurate energy results for a-bonded ligands, with less than 10% error 
in the ligand-field stabilization energy (defined below). 

This approach is extended in the present paper, in two ways. First, we develop energy 
functions which treat 7r-bonded ligands as well. This results in new and distinct angular 
forms. Second, we present results for the lowest-energy structures of small clusters, which 
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are a guide to understanding the bonding preferences of transition metals in proteins. 

The organization of the remainder of the paper is as follows. Section II develops the 
formalism underlying the angular forces. Section III gives tests of the functional form of the 
angular forces by comparing with results from diagonalization of simple cluster Hamiltonians. 
Section IV gives the small-cluster results for ideal examples, and for an analogue of the copper 
environment in blue-copper proteins. 



II. MODEL 

We model the electronic structure of the environment of a single transition metal ion 
in terms of orbitals \L,/i) localized on the ligands and \M, v) localized on the single metal 
ion. The orbitals taken to be orthogonal for simplicity of calculation. Here "L" denotes a 
particular ligand atom and the \L,fj) are distinct orbitals on that atom. The Hamiltonian 
takes the following form: 

H = E £ lA L i ( L > A*l + E £ m,u\M, v) (M, u\ 

+ Y,[KM\L,v){M,v\ + h 1 Z} L \M,is){L, f ,\] . (1) 

L,n,v 

Here the e terms are on-site energies for the orbitals, and the h terms are hybridization 
energies between the ligand orbitals and the metal orbitals. We ignore electron-interaction 
terms as well as explicit electrostatic effects. The energy associated with this Hamiltonian 
can then be obtained by diagonalizing its matrix as given in the \L,fj) — \M, v) basis, and 
taking the sum of the energies of the occupied eigenvectors. 

To obtain a real-space description of the energetics of this Hamiltonian, we make the 
simplifying approximation that the dimensionless ratios \h^ M \j '{sm,v — £ l,ix) are small and 
can be used as expansion parameters. This will be the case if the bonding is primarily ionic 
as opposed to covalent. In the case of transition metal ions, this allows a great simplification 
by the application of ligand-field theory, which gives an effective Hamiltonian for the ci-shell: 

H d = Y,h7W,n)(MM , (2) 

where 

K" = E( £d ~~ £ L,r,)~ 1 h l MLhTM ■ ( 3 ) 

L,T) 

Note that we have assumed all of the d-orbitals on the transition metal ion to have the same 
energies before ligand-field effects are "turned on" . 

In a fully quantum-mechanical ligand-field-theory calculation, one would numerically 
diagonalize the matrix Hd- In order to obtain a force field with a nearly classical form, we 
instead work with the moments of this matrix, defined in terms the traces of its powers, as 
follows. The first moment is defined as the average energy of the rf-complex, or 

e = (l/5)M d = (1/5) E(^ - e^)' 1 E MMV 
= (l/5)E( e d- e O"Vro ( 4 ) 
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where 

A* 

Thus e is given as a sum of independent contributions from the ligands. Because the ci-shell 
by itself has spherical symmetry, the contribution from each ligand is a radial function (no 
angular dependence) of the metal-ligand distance. Our major interest is in the angular terms 
resulting from the ligand-field splitting, so we do not consider the e term further. 

The width W of the (i-complex corresponding to the ligand-field splitting can in the 
simplest picture be described in terms of the second moment or variance 5s 2 of the eigen- 
values of Hd- We expect that W oc \J 5s 2 . We can obtain a simple real-space form for 5s 2 as 
follows: First we note that 5s 2 = (1/5) Y, n { £ n — £ ) 2 = (l/5)Tr(^ — el) 2 where the e n are 
the eigenvalues of H d . From Eqs. (|2]), (|3|), and (|5|), we see that 

VH 2 = £ h?h? 

= E i £ d £ L,r)) ^( £ d £ L',r)')~^ ^L'M^'ML ^LM^ML' 

L,n,L',n' n v 

— E { £ d — £ L,rj) { £ d — £ L',r)>)~ 9%f ■ (6) 
L,r),L',r)' 

Then the variance is 

5e 2 = (l/5)T[(H d - Elf = (1/5) (m, 2 - 5^ 2 ) 

= (1/5) J2 ( £ d- £ L, v )~ 1 ( £ d- £ L', v >)~ 1 [g 2 . ri > ~ (V%w^w] • (7) 

L,n,L',rj' 

We will focus on the ligand-field stabilization energy, which is defined as the sum of the 
energies of the occupied orbitals relative to s, or E LFSE = J2' n ( £ n — £ ), where the sum is over 
only the occupied eigenfunctions. We assume that -Exfse is proportional to the ligand-field 
splitting W , and also identify W with the standard deviation V 5s 2 . Thus our expression 
for the ligand-field stabilization energy has the form 

r _ r , -I 1/2 

#LFSE = ~\ a/5 ( £ d ~ £ L,n) ( £ d - £ L', V ') [9^' ~ (V^SW&jwJ f > ( 8 ) 

^ L,r),L',ri' ' ' 

where a is a dimensionless constant which will later be used as a fitting parameter. 

The accuracy of this form will be demonstrated in the next section via specific numerical 
experiments. At present, we will show how this form for -Elfse results in an expression for 
the energy in terms of simple angular interactions between he ligands. We thus need to 
develop analytic forms for the g m >. First we note that we can write 

gnr f = {L,ri\HPlH\L\Tf) 

= W> , (9) 

where Pd = \M, fj)(M, fj,\ is the projection onto the (i-subspace of the transition metal 
ion Here we define |\&) = PdH\L,r]) and we have have used the relation Pj = Pd, which 
follows from P^'s being a projection operator. 



4 



Consider first the case where both the \L,rj) and \L f ,i] f ) orbitals have a-character with 
respect to the metal ion. Then they couple only to the ci-orbitals \M,a) and \M,a') that 
have ex-character with respect to the bond axes, so that \^/) = h a \M, a) and |\&') = h' a \M, a'), 
where h a and h' a are the appropriate coupling strengths. Thus g m ' = ho-h^M, a\M, a'), and 
the angular dependence is contained in the last inner product. But this is simply a matrix 
element of a rotation about M which carries L into L'. We note that \M, a) is equivalent to 
\M,m = 0), where m = denotes the angular dependence of the spherical harmonic Y^m- 
Choosing our coordinate system so that L is along the z-axis and L' is in the z-x plane, we 
find that 

(M,a\M,a') = D$(0,e,0) 

= (3cos 2 #- l)/2 (10) 

where the D-term is a matrix element of the / = 2 representation of the rotation group, and 
the second equality follows from the explicit formulas of the D-terms given in Ref. (Wigner, 
1959). In summary, for cr-bonded ligands, 

g m < =/iX(3cos 2 #-l)/2 . (11) 

For the case of a 7r-bonded orbital L and a <T-bonded orbital L', we consider only the case 
in which the axis of the orbital lies along the circle connecting L and L'; if it is perpendicular 
to this circle, L and L' have different inversion symmetries so their coupling vanishes. We 
write |#) = K\M,ir} where \M,n) = (l/ v / 2)(-|M,m = l) + \M,m = -1)) and m the usual 
aximuthal angular momentum index for the spherical harmonics. (The "— " sign comes 
from the definition of the spherical harmonics). Then, again choosing a coordinate system 
in which L is along the z-axis and V is in the z-x plane, and following reasoning parallel to 
the a-bonded case, we see that 

(M,n\M,a') = (l/y/2) [-£>{? (0, 0, 0) + D%(0, 9, 0) 

= v / 3sin^cos^ . (12) 

and 

g m < = Vsh^h^ sin 6 cos 6 . (13) 

Finally we turn to the case of two 7r-bonded orbitals. To obtain the subsequent results, 
it is sufficient to consider the case 1) in which the orbitals are parallel to the arc connecting 
L and L f , and the case 2) in which they are perpendicular to it. By reasoning similar to 
that above, one sees that in case 1) 

(M, n\M, tt') = (1/2) [D$(0, 9, 0) + D^_i(0, 9, 0) - D?\(0, 9, 0) - £^(0, 9, 0) 

= cos2^ (14) 

and 

g m < = KK cos 261 • ( 15 ) 
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In case 2), 



(M, ir\M, tt') = (1/2) [D$(0, 9, 0) + (0, 0, 0) + £>fi(0, 0, 0)/ + dL^O, 0, 0)" 

= cos# (16) 

and 

g m < = h n h' n cos 9 . (17) 

Thus, combining Eqs. (§), (Q), (HH), dT5j) , and (|17D , we find the following explicit form 
for .Elfse as a sum of ligand-ligand interactions: 



£lfse = -/ E U(L, m L', V ') (18) 

where the ligand-ligand interaction U is given as follows: 
For a-a interactions, 

U(L, V ; L', r/) = (a/5)e a (r)e a (r')u aa (9) (19) 

for tx-u interactions, 

U(L, m L', i) = (a/5)e n (r)e a (r')u^(9) (20) 

and for ir-ir interactions, 

U(L, v, L', rf) = (a/h)e^r)e^r')u^(9) (21) 

where 

e CT (r) = K 2 /(£ d -e LjV ), (22) 
eir(r) = K 2 /(e d - e L rf), (23) 
Uaa {B) = (9 cos 4 9-6 cos 2 9 + l/5)/4, (24) 
u na (9) = (-3 cos 4 9 + 3 cos 2 9 - 2/5), (25) 
u^{9) = (4cos 4 ^-3cos 2 ^ + l/5) . (26) 

In each case, C/ is given as a product of radial terms involving the two ligands, and a simple 
angular function. (Note that although denoted an interaction here, U does not have units 
of energy because of the square root in Eq. (|IBD). In the calculation of the terms involving 
Ti-cx interactions, each term involves a sum over two 7r-orbitals, parallel and perpendicular 
to the arc connecting the two ligands. In the latter case, g m > in Eq. (H) vanishes, but g m 
and g v ' v ' do not. In the calculation of the terms involving it-it interactions, one has a similar 
scenario except that one sums over two pairs of 7r-orbitals. 

These forms are plotted out in Figure la. Note that the a-a interactions are fairly simi- 
lar in form to the 7r-7r interactions, both having pronounced minima at 180° (as well as the 
physically irrelevant one at 0°), and a shallower minimum at 90°. The <j-tt interaction is 
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complementary to these, having minima at 45° and 135°. We shall see later that these differ- 
ences lead to large differences in ground-state structures of small clusters. For comparison, 
we show in Figure lb an empirical angular interaction curve (Comba et al., 1995) assumed 
in some previous calculations of small transition-metal structures. The angular dependence 
is based on the observed square structure of small complexes of the transition metals of 
interest, and is proportional to sin 2 29. Since the metals which have square coordination 
generally have predominantly cr-bonds to their neighbors, the most relevant comparison is 
to the a-a curve in Figure la. We see that the behavior is quite different. The empiri- 
cal curve has equivalent minima at 180° and 90°, while in the theoretical curve the 180° 
minimum is much deeper. 
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III. TESTS OF FUNCTIONAL FORM 



In order to evaluate the accuracy of the semiclassical form of Eq. (H) for the ligand- 
field stabilization energy, we have performed explicit tests for small clusters. These clusters 
consist of a central transition-metal ion with four neighbors placed at random orientations 
at random distances relative to the central ion. We consider only the minority-spin orbitals, 
as these determine -Exfse for high-spin late transition metals. Four electrons are placed 
in these states, corresponding to Cu 2+ ; other band filling values gives similar results. The 
random distances are taken into account by varying the couplings h a and h v uniformly over 
a finite interval ranging from zero to /i^ ax or /i^ ax . We take /?^ ax = 0.5/i^ ax . The energies 
are obtained by explicit diagonalization of a tight-binding Hamiltonian for this cluster. All 
of the ligand orbitals are taken to have the same value of El. v - The only dimensionless 
variable that enters the results is then 7 = \h max /(ed — £l,ti)\- For small values of 7, the 
bonding is primarily ionic, and for larger values it acquires more covalent character. We use 
7 = 0.1 for our results, but even when 7 is significantly larger we find that the accuracy of 
the semiclassical form is essentially the same. 

We have fitted the energies of these clusters to the following semiclassical form for the 
energy: 

£lfse= E U(L,r,;L',Tf)+P , (27) 

L,V, L ',V' 

involving two parameters a and /3 (where U contains a). We use a database of 10,000 clusters 
to determine the parameters, and then test them on a set of 1000 clusters not included in the 
"training" set. Typical results are shown in Fig. 2, which shows the semiclassical energies 
vs. exact energies for 1000 clusters. Figure 2a corresponds to three a-bonded ligands and 
one 7r-bonded one, while Fig. 2b corresponds to four 7r-bonded ligands. The rms errors in 
Elfse for these two cases are 9% and 14%, respectively. For clusters with two and three 
7r-bonded ligands, the rms errors are 10% and 11%. This is to be compared with rms errors 
of 25%-30% that are obtained with empirical force fields (Carlsson, 1998). 
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IV. SMALL CLUSTER MINIMUM-ENERGY STRUCTURES 



In this section, we describe some of the implications of the angular forms developed 
above for the structure of small model clusters consisting of a metal atom and four ligands. 
Such cluster calculations cannot treat the protein environment accurately; this would await 
parametrization and incorporation of the force field into protein codes, which is in progress 
in our group. However, from the small-cluster calculations it is possible to see the structural 
preferences of the angular forces by themselves. These are an important factor in the final 
structure adopted by a protein. 

To keep the calculations as simple as possible, we place the ligands at frozen bond 
lengths from the central metal ion. Here "frozen means" that they are fixed at a given set 
of values, but these values are not necessarily the same for all of the ligands. The energy 
terms include the ligand-field stabilization energy as described by our angular terms, as well 
as a repulsive radial interaction between ligands forbidding close approaches. The latter 
term is needed because otherwise spurious structures involving 45° bond angles can appear 
when 7r-bonding ligands are present. For simplicity, we assume that the values of h a and 
h n at a given distance are the same. To evaluate the magnitudes of these couplings, we 
assume a value of 100 kJ/mole = 23.9 kcal/mole = 1.04 eV for the ligand-field stabilization 
energy of a cluster with four cr-ligands in square coordination; this number would correspond 
to the more strongly stabilized complexes (Cotton, 1972). For clusters with unequal bond 
lengths, we assume an exponential decay, so that h a (r) = h a (r ) exp [—k(t — r )], where r 
is the reference distance and k is a decay parameter. Since the most prominent case that 
we consider of a distant ligand is sulfur, we identify k with the spatial decay rate coming 
from the measured first ionization energy E 1 of sulfur, 10.4 eV (= 240 kcal/mole), using 
the formula h 2 K 2 /2m = E x . This yields k = 1.65A -1 . The ligand-ligand terms contain an 
exponential term taken as the repulsive part of the van der Waals interactions as given in the 
"MM2" force field (Sprague et al., 1987). In our model clusters, we use parameters and bond 
lengths typical for copper or nickel interacting with nitrogen and/or sulfur ligands. This 
is because copper and nickel have the largest ligand-field energies among the 3d transition 
metals, and typical ligands for these metals are nitrogen and sulfur. We have examined four 
simple geometries; the corresponding minimum-energy clusters are shown in Fig. 3. 

a) A cluster with all four ligands cr-bonded and placed at equal distances of 2. OA. We 
take repulsion parameters appropriate for nitrogen. The result is a square-planar cluster (cf. 
Fig. 3a), consistent with the observed structures of many copper and nickel complexes. With 
this parameterization, the square structure is quite strongly favored over the tetrahedral 
one, by 0.31 eV or 7.0 kcal/mole. The square structure remains stable under increases of 
the ligand-ligand repulsion strengths, of up to a factor of over three. We note, however, 
that the stability of the square structure is strongly sensitive to the choice of repulsive 
terms. For example, we find that when the van der Waals term from the OPLS-II force 
field (Jorgensen et al., 1988) is used instead of the MM2 form, an increase of the repulsion 
strengths of only 20% is enough to stabilize the tetrahedral structure. We expect that input 
from fully quantum-mechanical total-energy calculations will be necessary to precisely pin 
down the competition between the electronic energy and the repulsion terms in determining 
the structural energy differences. 

b) A cluster similar to that of a), but with one 7r-bonded ligand. This leads to a sub- 
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stantial deviation from planarity, as shown in Fig. 3b, where the 7r-bonded ligand is atom 4. 
As expected from the maximum of the n-a potential in Fig. 1, the angle between atom 4 
and atoms 1-3 is greater than 90°, about 123°. The angle between the u-bonded atoms 1, 

2, and 3, is slightly greater than the ideal value of 90°, because of the repulsive energy term. 

c) A cluster with with two cr-bonded nitrogen ligands 1 and 2, one a-bonded sulfur ligand 

3, and a 7r-bonded sulfur ligand 4. The distances are 2.04A for the nitrogen ligands, 2.18A 
for the 7T-bonded sulfur ligand, and 2.64A for the u-bonded sulfur ligand. This geometry 
is motivated by the observed geometry of copper sites in "blue-copper" proteins. In these 
proteins, the nitrogen ligands and a cysteine sulfur ligand are close to the copper, and it 
is believed that the cysteine sulfur ligand is primarily 7r-bonded. An additional cr-bonded 
methionine sulfur ligand is farther from the copper. The values of the ligand distances here 
are taken from cluster calculations for blue-protein models (Ryde et al., 1996). The van der 
Waals parameters are taken from the OPLS parameter set but we do not have a reliable 
procedure for determining the differences between the ligand-ion interactions of the nitrogen 
and the sulfur ligands. Since the greater distance of the sulfur at 2.18A from the copper will 
be compensated to some extent by the greater size of sulfur relative to nitrogen, we simply 
assume that the three close ligands have the same coupling strength to the central atom. 
We have varied the ratio of the sulfur to nitrogen coupling strengths by up to 30% in both 
directions and found changes of only a few degrees in the bond angles of the cluster. For the 
far ligand, we use the scaling procedure described above, which leads to a coupling that has 
50% of the strength of the close ligands. We have varied this coupling from 0% to 75% of the 
close-ligand value, and again found bond-angle changes of only a few degrees. The lowest- 
energy structure for this cluster is indicated in Fig. 3c. The geometry is trigonal, with bond 
angles of 125° between the cysteine sulfur and the nitrogens; the nitrogen-nitrogen bond 
angle is 91°. By comparison, the optimal values for protein models obtained by quantum- 
chemical calculations (Ryde et al., 1996) are 125° for the cysteine sulfur-nitrogen bond angle 
and 103° for the nitrogen-nitrogen bond angle. The results obtained by the present method 
are do of course not have quantitative accuracy, but the overall structure of the coordination 
shell is quite similar to that obtained in (Ryde et al., 1996), which is shown in Fig. 3d. The 
stability of this structure is not predetermined by the bond lengths that we used as input. 
To demonstrate this, we have considered a cluster with the same bond lengths as the blue- 
protein fragment, with with all ligands a-bonded. In this case, the structure is planar, as 
shown in Fig. 3e. Thus the formation of the trigonal structure is directly related to the 
special character of the angular interactions associated with 7r-bonding. 
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V. CONCLUSION 



We have seen that it is possible to extract the functional form of angular forces around 
transition metals in biomolecules by using an approximate treatment of the ligand-field 
Hamiltonian. The approximations made here, involving primarily ionic bonding, are com- 
plementary to those used in (Landis et al., 1998). The resulting angular forms have simple 
trigonometric forms, which are very different for 7r-bonding as compared to cr-bonding. 
Ligand-field energies are represented well by these angular forces, and much better than by 
force fields with assumed functional forms. Future work should aim to include the functional 
forms derived here in widely used biomolecular simulation packages. 
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FIGURES 



FIG. 1. (a) Interaction between ligands of central transition metal ion. The dimensionless 
quantity u (cf. Eqs. (19)-(21)) is either w CT)CT (solid line), u W)CT (dashed line), or u n<7r (dotted line). 
The quantity — u is plotted so that negative values will correspond to lower energies, (b) Empirical 
force field from (Comba et al., 1995). 

FIG. 2. Energies obtained by semiclassical force field vs. exact energies for small model clusters. 
Energies given in units of h max , maximal coupling strength between ligands and transition metal, 
(a) One 7r-bonded ligand and three cr-bonded ligands. (b) Four 7r-bonded ligands. 

FIG. 3. Lowest-energy structures obtained in relaxations of five-atom cluster using angular 
forces, (a) Four equivalent cr-bonded neighbors, (b) Three cr-bonded neighbors and one 7r-bonded 
neighbor 4, all at the same distance, (c) Two cr-bonded nitrogen neighbors 1 and 2, a near 7r-bonded 
sulfur neighbor 3, and a far cr-bonded neighbor 4. (d) Structure for blue protein model obtained 
by quantum-chemical calculations (Ryde et al., 1996). (e) Three near cr-bonded neighbors and a 
far cr-bonded neighbor. 
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